# Analysis of Streamflow data
# Created as part of Contract GSA-IS17JHQ0177
# Author  Carl James Schwarz
#         StatMathComp Consulting by Schwarz, Inc.
#         cschwarz.statsfuca@gmail.com
# Change Log
#     2017-01-30 CJS First Edition

options(width=200)  # how wide to print output

# load the functions used in the analysis
# This will also load the relevant libraries -- see the loadlibraries file in the Functions folder for details
dummy <- lapply(file.path("../../Functions",list.files(path=file.path("../../Functions"), pattern = "(?i)[.]r$", recursive = TRUE)), 
       function(x){
         cat("Loading ",x,"\n");
         source(x)})
## Loading  ../../Functions/compare.annual.stat.R 
## Loading  ../../Functions/compare.frequency.with.hec.R 
## Loading  ../../Functions/compare.longterm.stat.R 
## Loading  ../../Functions/compare.percentile.longterm.stat.R 
## Loading  ../../Functions/compute.Q.percentile.longterm.R 
## Loading  ../../Functions/compute.Q.stat.annual.R 
## Loading  ../../Functions/compute.Q.stat.longterm.R 
## Loading  ../../Functions/compute.volume.frequency.analysis.R 
## Loading  ../../Functions/load_packages.r
## Warning: package 'car' was built under R version 3.3.2
## Loading required package: MASS
## Loading required package: survival
## Warning: package 'ggplot2' was built under R version 3.3.2
## Warning: package 'zoo' was built under R version 3.3.2
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
# Define variables that provide options for the analysis

Station.Code <- '08HA011'
Station.Name <- 'Upper Cowichan River'
Station.Area <- 826    # square km's

start.year   <- 1965  # when do you want the analysis to start at.
end.year     <- 2012  # what is last year with complete data

# Get the data
flow <- read.csv(file.path("08HA011_STREAMFLOW_SUMMARY.csv"), header=TRUE, as.is=TRUE, strip.white=TRUE)

# Create a date variable 
flow$Date <- as.Date(paste(flow$Year,flow$Month, flow$Day, sep="-"), "%Y-%m-%d")

# basic structure of flow
str(flow)
## 'data.frame':    42734 obs. of  7 variables:
##  $ Day  : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Month: int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Year : int  1900 1900 1900 1900 1900 1900 1900 1900 1900 1900 ...
##  $ DOY  : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Q    : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Hour : chr  "24:00" "24:00" "24:00" "24:00" ...
##  $ Date : Date, format: "1900-01-01" "1900-01-02" "1900-01-03" "1900-01-04" ...
# Some preliminary screening

# Are there any illegal dates that could not be created?
cat("Number of illegal date values ", sum(is.na(flow$Date)), "\n")
## Number of illegal date values  0
flow[ is.na(flow$Date),]     # sorry Microsoft, but feb 29 1900 is NOT a leap year
## [1] Day   Month Year  DOY   Q     Hour  Date 
## <0 rows> (or 0-length row.names)
dim(flow)  # size before removing illegal dates
## [1] 42734     7
flow <- flow[ !is.na(flow$Date),]
dim(flow)  # size after removing illegal dates
## [1] 42734     7
# Table of simple statistics by year to help check for outliers, etc
flow.sum <- plyr::ddply(flow[ flow$Year >= start.year & flow$Year <=end.year,], "Year", plyr::summarize,
         n.days   = length(Year),
         n.Q      = sum (!is.na(Q)),
         n.miss.Q = sum ( is.na(Q)),
         min.Q    = min (Q, na.rm=TRUE),
         max.Q    = max (Q, na.rm=TRUE),
         mean.Q   = mean(Q,na.rm=TRUE),
         sd.Q     = sd  (Q,na.rm=TRUE))
flow.sum
##    Year n.days n.Q n.miss.Q min.Q max.Q   mean.Q     sd.Q
## 1  1965    365 365        0  3.54 196.0 47.72625 44.55059
## 2  1966    365 365        0  5.80 362.0 60.32545 65.58664
## 3  1967    365 365        0  4.11 213.0 61.70992 53.65439
## 4  1968    366 366        0  5.10 450.0 69.52503 71.15098
## 5  1969    365 365        0  4.98 150.0 47.50655 35.52722
## 6  1970    365 365        0  4.25 153.0 36.84767 34.33385
## 7  1971    365 365        0  6.00 187.0 55.62266 39.62329
## 8  1972    366 366        0  5.83 425.0 53.93230 59.41003
## 9  1973    365 365        0  3.11 337.0 53.33545 64.81365
## 10 1974    365 365        0  3.28 343.0 69.46066 64.42318
## 11 1975    365 365        0  6.43 345.0 66.59186 65.38940
## 12 1976    366 366        0  4.62 165.0 47.77153 37.67202
## 13 1977    365 365        0  4.16 233.0 47.69573 47.43126
## 14 1978    365 365        0  4.64 101.0 30.25252 21.66569
## 15 1979    365 365        0  3.84 336.0 43.46542 56.16374
## 16 1980    366 366        0  4.42 385.0 58.65967 62.00653
## 17 1981    365 365        0  4.66 244.0 60.31227 54.73371
## 18 1982    365 365        0  5.32 217.0 49.76348 46.18497
## 19 1983    365 365        0  5.71 303.0 61.29674 69.73937
## 20 1984    366 366        0  4.12 197.0 58.12112 41.60186
## 21 1985    365 365        0  3.36  86.4 27.75455 20.36948
## 22 1986    365 365        0  3.45 400.0 52.85359 63.78841
## 23 1987    365 365        0  3.60 198.0 47.36142 50.49256
## 24 1988    366 366        0  4.14 139.0 46.11101 37.02558
## 25 1989    365 365        0  4.08 176.0 40.17173 35.15942
## 26 1990    365 365        0  4.89 326.0 59.71548 63.55819
## 27 1991    365 365        0  4.61 248.0 45.93523 47.93800
## 28 1992    366 366        0  3.74 324.0 42.74970 56.93322
## 29 1993    365 365        0  4.89 191.0 33.46789 33.87104
## 30 1994    365 365        0  4.80 257.0 47.67293 52.08854
## 31 1995    365 365        0  4.67 286.0 71.34422 77.38259
## 32 1996    366 366        0  4.63 254.0 50.60377 45.75725
## 33 1997    365 365        0  7.45 359.0 79.12233 55.64481
## 34 1998    365 365        0  3.96 325.0 57.35233 71.50327
## 35 1999    365 365        0  6.01 282.0 71.80063 57.92698
## 36 2000    366 366        0  5.23  88.7 36.60932 22.77473
## 37 2001    365 365        0  4.21 255.0 41.97326 42.51711
## 38 2002    365 365        0  4.98 289.0 49.30463 50.89087
## 39 2003    365 365        0  2.68 312.0 66.03233 67.18302
## 40 2004    366 366        0  2.77 213.0 43.52918 43.62734
## 41 2005    365 365        0  4.96 321.0 48.51677 53.08182
## 42 2006    365 365        0  2.58 333.0 70.73003 85.25151
## 43 2007    365 365        0  5.64 381.0 71.45077 67.47710
## 44 2008    366 366        0  5.30 166.0 40.58847 30.46670
## 45 2009    365 365        0  5.48 407.0 53.16282 68.14245
## 46 2010    365 365        0  4.39 365.0 64.36416 61.11077
## 47 2011    365 365        0  5.68 222.0 57.97247 47.16704
## 48 2012    366 366        0  3.55 229.0 58.51451 48.18720
# visuallize the min, max, and mean
plotdata <- reshape2::melt(flow.sum,
             id.var='Year',
             measure.var=c("min.Q","max.Q","mean.Q","sd.Q"),
             variable.name='Statistic',
             value.name='Value')
ggplot2::ggplot(data=plotdata, aes(x=Year, y=Value))+
  ggtitle("Summary statistics about Q over time")+
  theme(plot.title = element_text(hjust = 0.5))+
  geom_line()+
  facet_wrap(~Statistic, ncol=2, scales="free_y")

# Some preliminary plots to check for outliers etc
ggplot2::ggplot(data=flow[flow$Year >= start.year,], aes(x=Date, y=Q))+
   ggtitle("Q by date")+
   theme(plot.title = element_text(hjust = 0.5))+
   geom_line(aes(group=Year))

ggplot2::ggplot(data=flow[flow$Year >= start.year,], aes(x=Date, y=log(Q)))+
   ggtitle("log(Q) by date")+
   theme(plot.title = element_text(hjust = 0.5))+
   geom_line(aes(group=Year))

ggplot2::ggplot(data=flow[flow$Year >= start.year,], aes(x=Date, y=log(Q)))+
   ggtitle("Q by date")+
   theme(plot.title = element_text(hjust = 0.5))+
   geom_line()+
   facet_wrap(~Year, scales="free_x")
## Warning: Removed 1461 rows containing missing values (geom_path).

# Compute the statistics on an annual basis


na.rm=list(na.rm.global=TRUE)

# Create the directory for the results of the analysis
# If you are comparing to the spreadsheet, the spreadsheet must also be in this directory.
# Similarly, if you are comparing to the HEC-SSP results, the vfa.rpt file must also be in this directory
#report.dir <- file.path("08HA011","MinFlows") # directory for output files and save objects
report.dir <- '.'  # current director
dir.create(report.dir)
## Warning in dir.create(report.dir): '.' already exists
cat("Reports and saved files will be found in ", report.dir, "\n")
## Reports and saved files will be found in  .
stat.annual <- compute.Q.stat.annual(Station.code=Station.Code, 
                          Station.Area=Station.Area, 
                          flow=flow, 
                          start.year=start.year, 
                          end.year=end.year,
                          write.stat.csv=TRUE,        # write out statistics 
                          write.stat.trans.csv=TRUE,  # write out statistics in transposed format
                          report.dir=report.dir,
                          na.rm=na.rm)


names(stat.annual)
## [1] "Q.stat.annual"       "Q.stat.annual.trans" "dates.missing.flows" "file.stat.csv"       "file.stat.trans.csv" "na.rm"               "Version"             "Date"
head(stat.annual$Q.stat.annual)
##   Year SW_1_day_MIN SW_1_day_MINDOY SW_3_day_MIN SW_3_day_MINDOY SW_7_day_MIN SW_7_day_MINDOY SW_30_day_MIN SW_30_day_MINDOY SW_ANNUAL_MIN SW_ANNUAL_MAX SW_ANNUAL_MEAN SW_ANNUAL_TOTALQ
## 1 1965         3.54             276     3.746667             276     4.224286             276      5.070000              223          3.54           196       47.72625       1505094910
## 2 1966         5.80             215     5.936667             224     6.305714             225      6.445667              244          5.80           362       60.32545       1902423451
## 3 1967         4.11             243     4.110000             256     4.320000             269      4.652333              270          4.11           213       61.70992       1946083966
## 4 1968         5.10             219     5.223333             220     5.258571             223      6.569000              237          5.10           450       69.52503       2198548223
## 5 1969         4.98             228     4.980000             238     5.071429             240      5.357667              256          4.98           150       47.50655       1498166495
## 6 1970         4.25             203     4.343333             204     4.642857             206      5.035667              230          4.25           153       36.84767       1162028161
##   SW_ANNUAL_YIELDMM SW_JFM_TOTALQ SW_AMJ_TOTALQ SW_JAS_TOTALQ SW_OND_TOTALQ SW_JFM_YIELDMM SW_AMJ_YIELDMM SW_JAS_YIELDMM SW_OND_YIELDMM JAN_MIN_DAILY_SW FEB_MIN_DAILY_SW MAR_MIN_DAILY_SW
## 1          1823.397     565669441     232522271      48689856     658213342       684.8298       281.5040       58.94656       796.8685             28.9             96.3             31.7
## 2          2304.754     732991673     314576353      74649600     780205824       887.3991       380.8430       90.37482       944.5591             57.2             54.4             46.7
## 3          2357.648     874480320     258498431      44503776     768601439      1058.6929       312.9521       53.87866       930.5102            106.0             58.6             49.6
## 4          2656.226    1148152321     227285567      93315456     729794880      1390.0149       275.1641      112.97271       883.5289             67.7             90.0             74.8
## 5          1815.003     443759040     540264386      76129631     438013438       537.2385       654.0731       92.16662       530.2826             30.0             28.3             36.5
## 6          1407.777     492134399     222619105      46146240     401128417       595.8044       269.5147       55.86712       485.6276             39.1             48.1             39.9
##   APR_MIN_DAILY_SW MAY_MIN_DAILY_SW JUN_MIN_DAILY_SW JUL_MIN_DAILY_SW AUG_MIN_DAILY_SW SEP_MIN_DAILY_SW OCT_MIN_DAILY_SW NOV_MIN_DAILY_SW DEC_MIN_DAILY_SW JAN_MAX_DAILY_SW FEB_MAX_DAILY_SW
## 1             28.9             8.07             6.51             4.81             4.50             4.11             3.54             72.8             68.8             92.6            196.0
## 2             43.0             7.36             6.77             6.09             5.80             6.26             6.37             30.6            119.0            260.0            119.0
## 3             34.5            28.10             6.23             5.80             4.11             4.11             4.70             51.0             49.3            213.0            187.0
## 4             36.8            10.10             8.07             6.40             5.10            12.60            14.40             72.8             88.1            450.0            211.0
## 5             81.3            58.00             6.51             6.09             4.98             5.24            24.40             23.5             34.0             98.5             65.7
## 6             24.7            12.50             5.47             4.25             4.70             5.24             6.23             19.5             45.0            150.0             80.7
##   MAR_MAX_DAILY_SW APR_MAX_DAILY_SW MAY_MAX_DAILY_SW JUN_MAX_DAILY_SW JUL_MAX_DAILY_SW AUG_MAX_DAILY_SW SEP_MAX_DAILY_SW OCT_MAX_DAILY_SW NOV_MAX_DAILY_SW DEC_MAX_DAILY_SW JAN_MEAN_SW FEB_MEAN_SW
## 1             86.7             70.8             56.1             15.9             9.23            11.40             9.66             67.1            171.0              188    45.10323   124.08571
## 2            117.0            140.0             43.9             31.7            26.90            13.30             8.75             43.3            182.0              362   123.84194    83.40357
## 3            138.0             77.6             39.1             32.6             7.39             6.74             6.51            167.0            156.0              187   143.38710   111.24643
## 4            135.0             76.5             39.6             20.8            16.30            16.00            19.30            122.0            141.0              177   204.73871   136.61724
## 5            103.0            141.0             79.0             64.8            12.70             6.09            36.20             39.1             49.6              150    60.32903    45.78214
## 6             72.2            138.0             23.8             11.8             5.66             5.32            11.60             27.2             71.4              153    71.65161    66.56071
##   MAR_MEAN_SW APR_MEAN_SW MAY_MEAN_SW JUN_MEAN_SW JUL_MEAN_SW AUG_MEAN_SW SEP_MEAN_SW OCT_MEAN_SW NOV_MEAN_SW DEC_MEAN_SW JAN_P50_SW FEB_P50_SW MAR_P50_SW APR_P50_SW MAY_P50_SW JUN_P50_SW JUL_P50_SW
## 1    54.01613    41.90333    38.00129    8.536333    5.498065    6.489677    6.397333    33.44871    96.13000   119.27097       39.6      117.0       54.9      32.70       40.5      7.930       5.15
## 2    74.49355    73.04000    31.00968   16.281000   13.913871    6.799032    7.396667    16.87290    61.63667   214.77419      111.0       81.7       69.7      72.80       35.7     11.000      14.00
## 3    82.62581    49.75333    32.96452   15.912667    6.551613    5.538710    4.676333    70.17258    94.94667   124.90645      124.0      103.0       81.0      46.70       33.1     11.295       6.51
## 4    96.12903    53.06333    21.15484   12.764000   11.491290    7.942258   15.920000    52.24194   100.08000   123.38065      203.0      124.0       98.3      50.15       19.0     11.800      10.90
## 5    64.00000   104.38000    69.13548   32.615333    8.174516    5.536774   15.202666    31.73226    40.38667    92.71935       61.2       42.6       73.1     103.00       70.5     32.850       7.65
## 6    51.97097    60.76333    17.22903    7.320333    5.243871    5.068064    7.147667    13.30613    43.72000    94.14839       62.0       65.8       46.7      55.35       18.7      6.555       5.32
##   AUG_P50_SW SEP_P50_SW OCT_P50_SW NOV_P50_SW DEC_P50_SW JAN_P80_SW FEB_P80_SW MAR_P80_SW APR_P80_SW MAY_P80_SW JUN_P80_SW JUL_P80_SW AUG_P80_SW SEP_P80_SW OCT_P80_SW NOV_P80_SW DEC_P80_SW JAN_P90_SW
## 1       5.24      5.915      25.70      95.10      111.0       33.7     101.40       38.2      29.70       34.3      7.360       4.96       4.96      4.960      19.40      79.90       89.2       29.7
## 2       6.54      7.335       7.82      54.10      182.0       91.2      59.70       54.4      51.70       23.4      6.928       6.48       6.03      6.698       6.94      35.76      165.0       66.8
## 3       5.38      4.450      73.10      98.55      118.0      118.0      77.62       64.0      39.10       30.3      7.220       6.29       5.10      4.296      47.90      65.54       86.7      113.0
## 4       6.48     15.350      38.20      89.45      121.0       90.9     106.00       79.0      39.32       13.0      8.964       8.78       5.47     13.460      21.20      80.24      106.0       76.2
## 5       5.66      7.180      32.80      42.20      104.0       36.8      34.94       40.5      89.68       64.6      6.880       6.65       5.10      5.320      25.20      36.44       40.2       33.7
## 6       5.18      5.540       9.23      48.85       88.1       47.0      60.60       40.5      37.22       12.9      6.190       4.98       4.70      5.240       7.31      22.72       62.3       44.7
##   FEB_P90_SW MAR_P90_SW APR_P90_SW MAY_P90_SW JUN_P90_SW JUL_P90_SW AUG_P90_SW SEP_P90_SW OCT_P90_SW NOV_P90_SW DEC_P90_SW SW_Oct_to_Sept_TOTALQ SW_Oct_to_Sept_YIELDMM SW_ONDJFM_TOTALQ
## 1      98.68       34.3      29.40       21.6      6.940       4.96       4.96      4.779       6.23      76.06       82.7            1132277554               1370.796       1143909313
## 2      58.00       51.3      47.37       13.7      6.770       6.48       6.03      6.577       6.65      31.64      162.0            1780430969               2155.485       1391205016
## 3      68.08       57.2      36.74       29.7      6.850       6.09       4.53      4.182       7.56      54.23       85.0            1957688352               2370.083       1654686144
## 4      95.08       74.8      38.20       12.8      8.580       7.67       5.10     12.690      16.10      78.70       96.0            2237354782               2708.662       1916753760
## 5      30.95       40.5      85.50       59.7      6.645       6.51       4.98      5.240      25.20      28.90       37.4            1789947937               2167.007       1173553920
## 6      57.29       39.9      33.26       12.6      5.993       4.53       4.70      5.240       6.60      20.47       48.4            1198913181               1451.469        930147837
##   SW_AMJJAS_TOTALQ SW_ONDJFM_YIELDMM SW_AMJJAS_YIELDMM
## 1        281212127          1384.878          340.4505
## 2        389225953          1684.268          471.2179
## 3        303002207          2003.252          366.8308
## 4        320601022          2320.525          388.1368
## 5        616394017          1420.767          746.2397
## 6        268765345          1126.087          325.3818
tail(stat.annual$Q.stat.annual)
##    Year SW_1_day_MIN SW_1_day_MINDOY SW_3_day_MIN SW_3_day_MINDOY SW_7_day_MIN SW_7_day_MINDOY SW_30_day_MIN SW_30_day_MINDOY SW_ANNUAL_MIN SW_ANNUAL_MAX SW_ANNUAL_MEAN SW_ANNUAL_TOTALQ
## 43 2007         5.64             197     5.680000             198     5.810000             245      5.952667              259          5.64           381       71.45077       2253271393
## 44 2008         5.30             218     5.330000             220     5.448571             221      5.669000              232          5.30           166       40.58847       1283504835
## 45 2009         5.48             242     5.533333             264     5.578571             268      5.625000              271          5.48           407       53.16282       1676542751
## 46 2010         4.39             242     4.460000             260     4.561429             242      4.792000              260          4.39           365       64.36416       2029788287
## 47 2011         5.68             256     5.840000             258     5.892857             260      6.210333              264          5.68           222       57.97247       1828219678
## 48 2012         3.55             281     3.596667             282     4.138572             277      4.687667              283          3.55           229       58.51451       1850369184
##    SW_ANNUAL_YIELDMM SW_JFM_TOTALQ SW_AMJ_TOTALQ SW_JAS_TOTALQ SW_OND_TOTALQ SW_JFM_YIELDMM SW_AMJ_YIELDMM SW_JAS_YIELDMM SW_OND_YIELDMM JAN_MIN_DAILY_SW FEB_MIN_DAILY_SW MAR_MIN_DAILY_SW
## 43          2729.800    1045543680     310594176      61282656     835850881      1265.7914       376.0220       74.19208      1011.9260             95.1             69.4             67.6
## 44          1550.696     555906242     292714561      77844672     357039360       673.0100       354.3760       94.24294       432.2510             53.5             42.8             46.8
## 45          2031.103     423973439     267626592      47885472     937057248       513.2850       324.0031       57.97273      1134.4519             32.8             27.3             28.5
## 46          2459.054     865581119     414793440      43291584     706122143      1047.9190       502.1712       52.41112       854.8694             99.6             64.7             53.3
## 47          2214.857     871732798     420056929      64668672     471761278      1055.3666       508.5435       78.29137       571.1396             91.7             65.1             58.1
## 48          2235.566     809023680     393163200      63979200     584203104       979.4476       475.9845       77.45666       707.2677             76.8             73.7             64.1
##    APR_MIN_DAILY_SW MAY_MIN_DAILY_SW JUN_MIN_DAILY_SW JUL_MIN_DAILY_SW AUG_MIN_DAILY_SW SEP_MIN_DAILY_SW OCT_MIN_DAILY_SW NOV_MIN_DAILY_SW DEC_MIN_DAILY_SW JAN_MAX_DAILY_SW FEB_MAX_DAILY_SW
## 43             55.0             25.0             6.53             5.64             5.65             5.77            21.50             50.8             71.8              304            107.0
## 44             28.1             19.5            15.00             5.59             5.30             5.80            16.70             27.5             30.8              166             65.5
## 45             38.6             16.6             7.03             5.92             5.48             5.50             5.50             63.9             83.2              121             49.2
## 46             70.0             17.6             7.82             5.14             4.39             4.39             4.41             54.2             62.3              365            113.0
## 47             50.9             29.9             6.38             6.26             6.09             5.68            12.90             35.5             37.9              211            130.0
## 48             68.9             22.2            15.10             6.44             5.53             4.08             3.55             57.1             63.4              229            133.0
##    MAR_MAX_DAILY_SW APR_MAX_DAILY_SW MAY_MAX_DAILY_SW JUN_MAX_DAILY_SW JUL_MAX_DAILY_SW AUG_MAX_DAILY_SW SEP_MAX_DAILY_SW OCT_MAX_DAILY_SW NOV_MAX_DAILY_SW DEC_MAX_DAILY_SW JAN_MEAN_SW FEB_MEAN_SW
## 43            278.0            105.0             67.4             24.5             9.42             9.06            19.20            112.0              157              381   176.93871    83.93214
## 44             71.0             47.2             83.4             49.9            15.70            14.20            18.20             32.4              111               56   101.85484    52.84483
## 45             90.2             68.0             63.3             16.4             7.03             6.50             5.76             84.2              407              239    73.64194    36.96071
## 46             99.3            133.0             66.7             67.4             7.89             5.40             9.38             85.8              108              261   176.43871    87.05357
## 47            222.0            107.0             74.5             50.4            11.50            10.40            26.80             51.7              193              135   135.17097    86.58214
## 48            132.0             94.1             85.3             25.4            17.60             7.12             5.99             68.6              130              194   128.37419    94.08966
##    MAR_MEAN_SW APR_MEAN_SW MAY_MEAN_SW JUN_MEAN_SW JUL_MEAN_SW AUG_MEAN_SW SEP_MEAN_SW OCT_MEAN_SW NOV_MEAN_SW DEC_MEAN_SW JAN_P50_SW FEB_P50_SW MAR_P50_SW APR_P50_SW MAY_P50_SW JUN_P50_SW JUL_P50_SW
## 43   137.61290    70.78333    33.91935    13.99467    6.841290    7.333226    8.996000    66.63226    94.61000   153.88065        149      83.40      146.0      70.10       29.9      14.50       6.12
## 44    56.26129    36.96333    50.99355    23.27333    9.100645    6.504194   13.907667    24.57742    70.35333    40.64194        104      51.20       54.8      37.85       56.2      21.10       7.74
## 45    51.26774    54.76333    36.80645    10.45433    6.449355    5.977419    5.633333    20.92806   208.24000   127.40645         76      35.80       41.3      57.10       33.2       9.27       6.39
## 46    68.10323    85.71667    43.28065    29.58833    5.936774    4.926452    5.476667    38.01000    75.67000   152.39677        163      84.55       64.5      78.95       43.5      20.50       5.65
## 47   112.09355    77.20667    53.54839    29.51900    8.702258    7.310323    8.403000    38.41935    66.76000    73.10968        128      80.50      100.0      81.10       55.7      30.90       8.92
## 48    85.66129    79.50000    50.99355    19.49000   12.847097    6.103548    5.101000    15.77774    89.49333   115.73226        122      90.40       85.1      74.85       48.4      18.00      13.90
##    AUG_P50_SW SEP_P50_SW OCT_P50_SW NOV_P50_SW DEC_P50_SW JAN_P80_SW FEB_P80_SW MAR_P80_SW APR_P80_SW MAY_P80_SW JUN_P80_SW JUL_P80_SW AUG_P80_SW SEP_P80_SW OCT_P80_SW NOV_P80_SW DEC_P80_SW
## 43       7.11      6.350       65.9      87.60      137.0      117.0      72.54       78.8      63.06       25.9      7.088       5.90       5.84      6.000      46.70      59.70      102.0
## 44       5.76     15.350       25.7      71.15       39.2       71.4      46.12       49.9      31.92       21.4     17.760       5.76       5.49     10.316      21.40      54.08       33.0
## 45       5.82      5.625       12.2     192.00      112.0       50.9      30.06       31.1      43.92       21.2      7.704       6.12       5.67      5.598       7.30     102.96       94.5
## 46       5.00      5.160       33.3      72.95      156.0      115.0      75.06       59.3      73.48       25.5      8.630       5.41       4.56      4.788       6.65      60.88       92.7
## 47       6.80      6.160       41.0      41.85       61.1      108.0      71.40       68.7      53.40       37.5     13.100       6.66       6.40      5.954      31.50      37.70       44.5
## 48       5.90      5.245       11.2      79.30       97.4       96.9      84.04       72.0      71.94       36.4     16.240      10.30       5.72      4.658       4.26      65.20       78.8
##    JAN_P90_SW FEB_P90_SW MAR_P90_SW APR_P90_SW MAY_P90_SW JUN_P90_SW JUL_P90_SW AUG_P90_SW SEP_P90_SW OCT_P90_SW NOV_P90_SW DEC_P90_SW SW_Oct_to_Sept_TOTALQ SW_Oct_to_Sept_YIELDMM SW_ONDJFM_TOTALQ
## 43      110.0      72.09       72.8      56.50       25.4      6.643       5.78       5.81      5.894      22.40      54.04       89.4            2401276320               2907.114       2029399488
## 44       58.8      45.54       48.5      29.67       20.7     15.840       5.69       5.36      6.635      20.50      34.40       31.9            1762316356               2133.555       1391757123
## 45       35.0      29.03       29.7      40.16       17.0      7.522       6.11       5.62      5.578       7.17      71.59       88.3            1096524863               1327.512        781012799
## 46      111.0      70.35       57.0      73.04       22.2      8.040       5.34       4.49      4.655       4.76      57.16       70.4            2260723392               2736.953       1802638367
## 47      103.0      67.00       63.6      52.03       34.3     10.061       6.42       6.27      5.895      31.40      37.23       41.0            2062580543               2497.071       1577854941
## 48       87.4      81.48       68.2      70.95       27.6     15.890       7.43       5.61      4.148       3.98      62.52       73.1            1737927359               2104.028       1280784958
##    SW_AMJJAS_TOTALQ SW_ONDJFM_YIELDMM SW_AMJJAS_YIELDMM
## 43        371876832         2456.9001          450.2141
## 44        370559234         1684.9360          448.6189
## 45        315512064          945.5361          381.9759
## 46        458085025         2182.3709          554.5824
## 47        484725601         1910.2360          586.8349
## 48        457142400         1550.5871          553.4412
head(stat.annual$Q.stat.annual.trans)
##                      Y1965      Y1966  Y1967      Y1968      Y1969      Y1970      Y1971      Y1972  Y1973      Y1974      Y1975      Y1976      Y1977  Y1978      Y1979      Y1980      Y1981
## SW_1_day_MIN      3.540000   5.800000   4.11   5.100000   4.980000   4.250000   6.000000   5.830000   3.11   3.280000   6.430000   4.620000   4.160000   4.64   3.840000   4.420000   4.660000
## SW_1_day_MINDOY 276.000000 215.000000 243.00 219.000000 228.000000 203.000000 223.000000 304.000000 285.00 308.000000 206.000000 223.000000 227.000000 212.00 187.000000 244.000000 252.000000
## SW_3_day_MIN      3.746667   5.936667   4.11   5.223333   4.980000   4.343333   6.293333   6.000000   4.47   3.443333   6.713333   4.976667   4.643333   4.70   3.920000   4.430000   4.833333
## SW_3_day_MINDOY 276.000000 224.000000 256.00 220.000000 238.000000 204.000000 224.000000 304.000000 253.00 309.000000 207.000000 223.000000 227.000000 213.00 188.000000 244.000000 249.000000
## SW_7_day_MIN      4.224286   6.305714   4.32   5.258571   5.071429   4.642857   6.718571   6.141429   4.57   3.712857   6.922857   5.405714   4.975714   4.78   3.941429   4.897143   4.931429
## SW_7_day_MINDOY 276.000000 225.000000 269.00 223.000000 240.000000 206.000000 225.000000 262.000000 279.00 309.000000 210.000000 225.000000 227.000000 215.00 177.000000 245.000000 253.000000
##                      Y1982      Y1983      Y1984      Y1985  Y1986      Y1987      Y1988      Y1989      Y1990      Y1991      Y1992  Y1993      Y1994      Y1995  Y1996  Y1997      Y1998      Y1999
## SW_1_day_MIN      5.320000   5.710000   4.120000   3.360000   3.45   3.600000   4.140000   4.080000   4.890000   4.610000   3.740000   4.89   4.800000   4.670000   4.63   7.45   3.960000   6.010000
## SW_1_day_MINDOY 213.000000 225.000000 203.000000 241.000000 297.00 221.000000 244.000000 280.000000 227.000000 233.000000 221.000000 276.00 233.000000 175.000000 239.00 231.00 264.000000 260.000000
## SW_3_day_MIN      5.360000   5.746667   4.420000   3.616667   3.49   3.733333   4.416667   4.120000   4.966667   4.706667   3.790000   5.06   4.880000   4.800000   4.65   7.52   3.996667   6.120000
## SW_3_day_MINDOY 234.000000 233.000000 203.000000 234.000000 297.00 222.000000 245.000000 281.000000 227.000000 210.000000 221.000000 247.00 234.000000 176.000000 239.00 236.00 249.000000 261.000000
## SW_7_day_MIN      5.381429   5.857143   5.021429   3.717143   3.68   3.968571   4.530000   4.368571   5.245714   4.765714   4.002857   5.09   5.037143   4.872857   4.72   7.59   4.018571   6.235714
## SW_7_day_MINDOY 219.000000 233.000000 207.000000 243.000000 298.00 314.000000 224.000000 284.000000 227.000000 214.000000 224.000000 251.00 245.000000 181.000000 243.00 236.00 250.000000 261.000000
##                      Y2000      Y2001      Y2002  Y2003      Y2004      Y2005      Y2006  Y2007      Y2008      Y2009      Y2010      Y2011      Y2012
## SW_1_day_MIN      5.230000   4.210000   4.980000   2.68   2.770000   4.960000   2.580000   5.64   5.300000   5.480000   4.390000   5.680000   3.550000
## SW_1_day_MINDOY 246.000000 230.000000 309.000000 271.00 228.000000 221.000000 256.000000 197.00 218.000000 242.000000 242.000000 256.000000 281.000000
## SW_3_day_MIN      5.300000   4.226667   5.106667   2.75   2.823333   5.020000   2.616667   5.68   5.330000   5.533333   4.460000   5.840000   3.596667
## SW_3_day_MINDOY 246.000000 232.000000 309.000000 277.00 229.000000 221.000000 257.000000 198.00 220.000000 264.000000 260.000000 258.000000 282.000000
## SW_7_day_MIN      5.391429   4.264286   5.317143   2.77   2.888571   5.131429   2.660000   5.81   5.448571   5.578571   4.561429   5.892857   4.138572
## SW_7_day_MINDOY 247.000000 232.000000 310.000000 277.00 233.000000 222.000000 259.000000 245.00 221.000000 268.000000 242.000000 260.000000 277.000000
head(stat.annual$dates.missing.flows)
## character(0)
stat.annual$file.stat.csv
## [1] "./08HA011-annual-summary-stat.csv"
stat.annual$file.stat.trans.csv
## [1] "./08HA011-annual-summary-stat-trans.csv"
# Compare the annual statistics with those in the Excel spreasheet
compare.annual <- compare.annual.stat(Q.filename=stat.annual$file.stat.csv,
                                      E.filename=file.path(report.dir,"08HA011_STREAMFLOW_SUMMARY.xlsx"),
                                      report.dir=report.dir,
                                      write.plots.pdf=TRUE,
                                      write.comparison.csv=TRUE)
## Warning: Removed 4 rows containing missing values (geom_point).
## Warning: Removed 4 rows containing missing values (geom_point).
names(compare.annual)
## [1] "stats.in.Q.not.in.E" "stats.in.E.not.in.Q" "diff.stat"           "plot.list"           "stat.not.plotted"    "file.comparison.csv" "file.plots.pdf"      "Version"            
## [9] "Date"
compare.annual$stats.in.Q.not.in.E
## character(0)
compare.annual$stats.in.E.not.in.Q
## [1] "SW_TOTAL_Q_CMS"    "Date_50P_Cummul_Q" "Date_25P_Cummul_Q" "Date_75P_Cummul_Q"
head(compare.annual$diff.stat)
##   Year Value.Q Value.E          diff mean        pdiff         stat
## 1 1965    3.54    3.54  1.469700e-10 3.54 4.151695e-11 SW_1_day_MIN
## 2 1966    5.80    5.80  2.651399e-10 5.80 4.571378e-11 SW_1_day_MIN
## 3 1967    4.11    4.11  4.856000e-10 4.11 1.181508e-10 SW_1_day_MIN
## 4 1968    5.10    5.10  3.674296e-10 5.10 7.204503e-11 SW_1_day_MIN
## 5 1969    4.98    4.98 -7.348966e-11 4.98 1.475696e-11 SW_1_day_MIN
## 6 1970    4.25    4.25  0.000000e+00 4.25 0.000000e+00 SW_1_day_MIN
names(compare.annual$plot.list)
## [1] "plot.daily"   "plot.yieldmm" "plot.mean"    "plot.per"     "plot.wyear"
l_ply(compare.annual$plot.list, function(x){plot(x)})

## Warning: Removed 4 rows containing missing values (geom_point).

## Warning: Removed 4 rows containing missing values (geom_point).

# Compute the long-term statistics

stat.longterm <- compute.Q.stat.longterm(Station.code=Station.Code, 
                          Station.Area=Station.Area, 
                          flow=flow, 
                          start.year=start.year, 
                          end.year=end.year,
                          write.stat.csv=TRUE,        # write out statistics 
                          write.stat.trans.csv=TRUE,
                          report.dir=report.dir)  # write out statistics in transposed format

names(stat.longterm)
## [1] "Q.stat.longterm"       "Q.stat.longterm.trans" "file.stat.csv"         "file.stat.trans.csv"   "na.rm"                 "Version"               "Date"
head(stat.longterm$Q.stat.longterm)
##   Month      mean median   max   min
## 1     1 110.46828  98.65 450.0 15.60
## 2     2  92.69174  82.90 303.0 14.70
## 3     3  77.38508  68.85 359.0 15.80
## 4     4  57.23731  52.40 154.0  8.67
## 5     5  35.27665  31.00 121.0  6.47
## 6     6  18.18226  13.60  73.4  2.99
tail(stat.longterm$Q.stat.longterm)
##    Month       mean  median   max   min
## 8      8   6.069792   5.550  87.9  2.77
## 9      9   8.989847   6.055  84.9  2.58
## 10    10  26.823743  14.200 312.0  2.72
## 11    11  85.113521  71.100 407.0  3.28
## 12    12 113.266129 104.000 425.0 24.20
## 13    13  53.262478  37.400 450.0  2.58
head(stat.longterm$Q.stat.longterm.trans)
##             JAN       FEB       MAR       APR       MAY      JUN       JUL       AUG       SEP       OCT       NOV      DEC  Longterm
## mean   110.4683  92.69174  77.38508  57.23731  35.27665 18.18226  9.739046  6.069792  8.989847  26.82374  85.11352 113.2661  53.26248
## median  98.6500  82.90000  68.85000  52.40000  31.00000 13.60000  6.755000  5.550000  6.055000  14.20000  71.10000 104.0000  37.40000
## max    450.0000 303.00000 359.00000 154.00000 121.00000 73.40000 64.800003 87.900002 84.900002 312.00000 407.00000 425.0000 450.00000
## min     15.6000  14.70000  15.80000   8.67000   6.47000  2.99000  2.850000  2.770000  2.580000   2.72000   3.28000  24.2000   2.58000
stat.longterm$file.stat.csv
## [1] "./08HA011-longterm-summary-stat.csv"
stat.longterm$file.stat.trans.csv
## [1] "./08HA011-longterm-summary-stat-trans.csv"
# Compare the longterm statistics with those in the Excel spreasheet
compare.longterm <- compare.longterm.stat(Q.filename=stat.longterm$file.stat.csv,
                                      E.filename=file.path(report.dir,"08HA011_STREAMFLOW_SUMMARY.xlsx"),
                                      report.dir=report.dir,
                                      write.comparison.csv=TRUE,
                                      write.plots.pdf=TRUE)
## Warning: Removed 2 rows containing missing values (geom_point).
names(compare.longterm)
## [1] "stats.in.Q.not.in.E" "stats.in.E.not.in.Q" "diff.stat"           "plot.list"           "stat.not.plotted"    "file.comparison.csv" "file.plots.pdf"      "Version"            
## [9] "Date"
compare.longterm$stats.in.Q.not.in.E
## character(0)
compare.longterm$stats.in.E.not.in.Q
## character(0)
head(compare.longterm$diff.stat)
##   Month   Value.Q   Value.E          diff      mean        pdiff stat
## 1     1 110.46828 110.46830 -2.042926e-05 110.46829 1.849333e-07 mean
## 2     2  92.69174  92.69174  3.771534e-07  92.69174 4.068900e-09 mean
## 3     3  77.38508  77.38508  6.253091e-07  77.38508 8.080487e-09 mean
## 4     4  57.23731  57.23731 -4.460229e-06  57.23731 7.792520e-08 mean
## 5     5  35.27665  35.27665 -3.460968e-06  35.27665 9.810931e-08 mean
## 6     6  18.18226  18.18226  3.881854e-06  18.18226 2.134968e-07 mean
names(compare.longterm$plot.list)
## [1] "plot.allstat"
l_ply(compare.longterm$plot.list, function(x){plot(x)})
## Warning: Removed 2 rows containing missing values (geom_point).

# Compute the long-term percentile statistics

percentile.longterm <- compute.Q.percentile.longterm(Station.code=Station.Code, 
                          Station.Area=Station.Area, 
                          flow=flow, 
                          start.year=start.year, 
                          end.year=end.year,
                          write.stat.csv=TRUE,        # write out statistics 
                          write.stat.trans.csv=TRUE,  # write out statistics in transposed format
                          report.dir=report.dir)

names(percentile.longterm)
## [1] "Q.percentile.stat"       "Q.percentile.stat.trans" "file.stat.csv"           "file.stat.trans.csv"     "na.rm"                   "Version"                 "Date"
head(percentile.longterm$Q.percentile.stat)
##   Month     P01    P02     P05    P10     P15    P20     P25    P30     P35    P40     P45   P50    P55   P60     P65   P70    P75    P80   P85   P90    P95     P99
## 1   Jan 321.000 282.26 239.650 194.00 169.900 152.60 141.000 129.00 120.550 112.00 106.000 98.65 92.215 84.94 78.2000 70.81 65.800 58.600 49.60 40.41 31.480 21.8740
## 2   Feb 245.450 229.80 186.000 154.00 133.000 122.00 114.000 107.00 101.000  95.40  89.150 82.90 77.275 73.10 68.5250 64.10 60.000 54.900 49.60 43.90 36.950 23.2550
## 3   Mar 209.520 188.26 161.650 132.30 112.950 101.00  92.900  85.40  79.000  74.84  71.385 68.85 65.115 62.60 58.9450 55.40 51.300 47.200 41.91 36.07 29.235 20.3960
## 4   Apr 131.610 120.00 107.000  93.91  85.000  77.34  72.425  68.00  63.700  59.70  56.000 52.40 48.655 45.66 43.6000 41.30 38.575 35.800 32.30 28.48 23.490 13.4000
## 5   May  84.082  79.13  72.800  65.83  60.295  54.30  49.000  43.70  40.055  36.62  33.400 31.00 29.000 26.40 23.7000 21.31 18.350 16.100 14.20 12.20  9.867  7.4122
## 6   Jun  63.954  57.20  49.005  40.40  32.215  26.82  23.025  20.63  18.800  16.84  15.200 13.60 12.000 10.70  9.0965  8.16  7.580  6.948  6.55  6.28  5.580  3.9500
tail(percentile.longterm$Q.percentile.stat)
##    Month     P01     P02      P05     P10     P15     P20    P25     P30      P35     P40     P45     P50    P55    P60     P65    P70     P75    P80     P85    P90     P95     P99
## 7    Jul  40.500  35.504  25.3650  17.500  14.395  12.300  10.80   9.799   8.5700   7.754   7.107   6.755  6.493  6.250  6.0790  5.850  5.6600  5.464  5.2500  4.960  4.3535  3.0900
## 8    Aug  13.313  11.700   9.6495   8.096   7.179   6.716   6.43   6.120   5.9400   5.800   5.670   5.550  5.460  5.368  5.2545  5.140  5.0000  4.860  4.6100  4.350  3.8800  3.0274
## 9    Sep  37.705  34.476  23.3150  16.700  14.200  11.920   9.35   8.263   7.3705   6.718   6.270   6.055  5.860  5.660  5.5000  5.340  5.2300  5.030  4.7900  4.390  3.8895  2.8639
## 10   Oct 157.130 125.000  90.9650  64.090  49.790  39.520  32.50  27.200  22.4000  19.000  15.985  14.200 12.800 11.580 10.2000  8.765  7.6175  6.674  5.7805  5.267  4.4235  3.1861
## 11   Nov 304.000 282.220 213.0500 161.100 134.150 120.000 111.00 103.000  95.5000  87.500  79.600  71.100 64.800 57.200 51.8950 46.900 42.5000 37.700 31.4000 20.090 10.7950  4.7238
## 12   Dec 306.520 269.520 231.0000 195.000 173.000 157.600 146.00 135.000 126.0000 118.000 111.000 104.000 94.900 88.100 81.3000 73.710 66.0750 58.600 51.4050 44.270 36.5000 26.7740
head(percentile.longterm$Q.percentile.stat.trans)
##        JAN    FEB    MAR    APR    MAY    JUN    JUL     AUG    SEP     OCT    NOV    DEC
## P01 321.00 245.45 209.52 131.61 84.082 63.954 40.500 13.3130 37.705 157.130 304.00 306.52
## P02 282.26 229.80 188.26 120.00 79.130 57.200 35.504 11.7000 34.476 125.000 282.22 269.52
## P05 239.65 186.00 161.65 107.00 72.800 49.005 25.365  9.6495 23.315  90.965 213.05 231.00
## P10 194.00 154.00 132.30  93.91 65.830 40.400 17.500  8.0960 16.700  64.090 161.10 195.00
## P15 169.90 133.00 112.95  85.00 60.295 32.215 14.395  7.1790 14.200  49.790 134.15 173.00
## P20 152.60 122.00 101.00  77.34 54.300 26.820 12.300  6.7160 11.920  39.520 120.00 157.60
percentile.longterm$file.stat.csv
## [1] "./08HA011-longterm-percentile-stat.csv"
percentile.longterm$file.stat.trans.csv
## [1] "./08HA011-longterm-percentile-stat-trans.csv"
# Compare the longterm percentile statistics with those in the Excel spreasheet
compare.percentile.longterm <- compare.percentile.longterm.stat(Q.filename=percentile.longterm$file.stat.csv,
                                      E.filename=file.path(report.dir,"08HA011_STREAMFLOW_SUMMARY.xlsx"),
                                      write.comparison.csv=TRUE,
                                      write.plots.pdf=TRUE,
                                      report.dir=report.dir)
names(compare.percentile.longterm)
## [1] "stats.in.Q.not.in.E" "stats.in.E.not.in.Q" "diff.stat"           "plot.list"           "stat.not.plotted"    "file.cmparsion.csv"  "file.plots.pdf"      "Version"            
## [9] "Date"
compare.percentile.longterm$stats.in.Q.not.in.E
## character(0)
compare.percentile.longterm$stats.in.E.not.in.Q
## character(0)
head(compare.percentile.longterm$diff.stat)
##   Month Value.Q Value.E          diff    mean        pdiff stat
## 1   Apr 131.610 131.610  1.136868e-13 131.610 8.638161e-16  P01
## 2   Aug  13.313  13.313 -4.482903e-10  13.313 3.367312e-11  P01
## 3   Dec 306.520 306.520  4.547474e-13 306.520 1.483581e-15  P01
## 4   Feb 245.450 245.450 -5.684342e-14 245.450 2.315886e-16  P01
## 5   Jan 321.000 321.000  0.000000e+00 321.000 0.000000e+00  P01
## 6   Jul  40.500  40.500  0.000000e+00  40.500 0.000000e+00  P01
names(compare.percentile.longterm$plot.list)
## [1] "plot.allstat"
l_ply(compare.percentile.longterm$plot.list, function(x){plot(x)})

# Compute the volume frequency analysis (similar to HEC SSP)

vfa.analysis <- compute.volume.frequency.analysis( Station.Code="08HA011", flow, 
                      start.year=1965, end.year=2016, use.water.year=FALSE, 
                      roll.avg.days=c(1,3,7,15),
                      use.log=FALSE,
                      use.max=FALSE,
                      fit.distr="PIII",
                      write.stat.csv=TRUE,
                      write.plotdata.csv=TRUE,
                      write.quantiles.csv=TRUE,
                      report.dir=report.dir)
## Warning in min(x$Q, na.rm = TRUE): no non-missing arguments to min; returning Inf

## Warning in min(x$Q, na.rm = TRUE): no non-missing arguments to min; returning Inf

## Warning in min(x$Q, na.rm = TRUE): no non-missing arguments to min; returning Inf

## Warning in min(x$Q, na.rm = TRUE): no non-missing arguments to min; returning Inf

## Warning in min(x$Q, na.rm = TRUE): no non-missing arguments to min; returning Inf

## Warning in min(x$Q, na.rm = TRUE): no non-missing arguments to min; returning Inf

## Warning in min(x$Q, na.rm = TRUE): no non-missing arguments to min; returning Inf

## Warning in min(x$Q, na.rm = TRUE): no non-missing arguments to min; returning Inf

## Warning in min(x$Q, na.rm = TRUE): no non-missing arguments to min; returning Inf

## Warning in min(x$Q, na.rm = TRUE): no non-missing arguments to min; returning Inf

## Warning in min(x$Q, na.rm = TRUE): no non-missing arguments to min; returning Inf

## Warning in min(x$Q, na.rm = TRUE): no non-missing arguments to min; returning Inf

## Warning in min(x$Q, na.rm = TRUE): no non-missing arguments to min; returning Inf

## Warning in min(x$Q, na.rm = TRUE): no non-missing arguments to min; returning Inf

## Warning in min(x$Q, na.rm = TRUE): no non-missing arguments to min; returning Inf

## Warning in min(x$Q, na.rm = TRUE): no non-missing arguments to min; returning Inf
names(vfa.analysis)
##  [1] "start.year"              "end.year"                "use.water.year"          "use.max"                 "roll.avg.days"           "Q.stat"                  "Q.stat.trans"           
##  [8] "plotdata"                "prob.plot.position"      "freqplot"                "fit.distr"               "fit"                     "fitted.quantiles"        "fitted.quantiles.trans" 
## [15] "file.stat.csv"           "file.stat.trans.csv"     "file.plotdata.csv"       "file.quantile.csv"       "file.quantile.trans.csv" "Version"                 "date"
vfa.analysis$file.stat.csv
## [1] "./08HA011-annual-vfa-stat.csv"
head(vfa.analysis$Q.stat)
##   Year  Measure    value
## 1 1965 Q001-avg 3.540000
## 2 1965 Q003-avg 3.746667
## 3 1965 Q007-avg 4.224286
## 4 1965 Q015-avg 4.633333
## 5 1966 Q001-avg 5.800000
## 6 1966 Q003-avg 5.936667
head(vfa.analysis$Q.stat.trans)
##   Year Q001-avg Q003-avg Q007-avg Q015-avg
## 1 1965     3.54 3.746667 4.224286 4.633333
## 2 1966     5.80 5.936667 6.305714 6.419333
## 3 1967     4.11 4.110000 4.320000 4.415333
## 4 1968     5.10 5.223333 5.258571 5.563333
## 5 1969     4.98 4.980000 5.071429 5.262000
## 6 1970     4.25 4.343333 4.642857 4.863333
vfa.analysis$freqplot

ggsave(plot=vfa.analysis$freqplot, 
       file=file.path(report.dir, paste(Station.Code,"-annual-vfa-frequency-plot.png",sep="")), h=6, w=6, units="in", dpi=300)

vfa.analysis$fitted.quantiles.trans
##    distr  prob Q001-avg Q003-avg Q007-avg Q015-avg
## 1   PIII 0.010 2.416806 2.581846 2.675951 2.838811
## 2   PIII 0.050 2.988016 3.145622 3.285016 3.450660
## 3   PIII 0.100 3.308550 3.462336 3.620229 3.787091
## 4   PIII 0.200 3.711903 3.861207 4.036056 4.204141
## 5   PIII 0.500 4.530930 4.672124 4.862202 5.031834
## 6   PIII 0.800 5.413735 5.547480 5.729219 5.899323
## 7   PIII 0.900 5.901163 6.031295 6.198934 6.368843
## 8   PIII 0.950 6.317356 6.444656 6.595474 6.764989
## 9   PIII 0.975 6.688462 6.813419 6.945789 7.114789
## 10  PIII 0.980 6.800741 6.925022 7.051201 7.220017
## 11  PIII 0.990 7.131889 7.254258 7.360602 7.528800
HEC.filename <- file.path(report.dir,"HEC-vfa-min.rpt")
compare.with.HEC <- compare.frequency.with.hec(Q.file.stat=vfa.analysis$file.stat.csv, 
                                               Q.file.plotdata=vfa.analysis$file.plotdata.csv,
                                               Q.file.quantile=vfa.analysis$file.quantile.csv,
                                               HEC.filename=HEC.filename,
                                               report.dir=report.dir,
                                               write.comparison.csv=TRUE,
                                               write.plots.pdf=TRUE)
names(compare.with.HEC)
## [1] "stats.in.Q.not.in.HEC" "stats.in.HEC.not.in.Q" "diff.stat"             "plot.list"             "stat.not.plotted"      "file.comparison.csv"   "file.plots.pdf"        "Version"              
## [9] "Date"
compare.with.HEC$stats.in.Q.not.in.HEC
## character(0)
compare.with.HEC$stats.in.HEC.not.in.Q
##  [1] "Q030-avg"   "Q060-avg"   "Q090-avg"   "Q120-avg"   "Q183-avg"   "Q030-pp"    "Q060-pp"    "Q090-pp"    "Q120-pp"    "Q183-pp"    "Q030-avg-q" "Q060-avg-q" "Q090-avg-q" "Q120-avg-q" "Q183-avg-q"
head(compare.with.HEC$diff.stat)
##   Year    Measure  value.Q value.HEC        diff     mean       pdiff
## 1 0.01 Q001-avg-q 2.416806     2.438 -0.02119401 2.427403 0.008731144
## 2 0.01 Q003-avg-q 2.581846     2.544  0.03784636 2.562923 0.014766873
## 3 0.01 Q007-avg-q 2.675951     2.612  0.06395081 2.643975 0.024187369
## 4 0.01 Q015-avg-q 2.838811     2.738  0.10081146 2.788406 0.036153798
## 5 0.05 Q001-avg-q 2.988016     2.976  0.01201559 2.982008 0.004029361
## 6 0.05 Q003-avg-q 3.145622     3.108  0.03762155 3.126811 0.012031925
names(compare.with.HEC$plot.list)
## [1] "plot.avg"   "plot.pp"    "plot.quant"
l_ply(compare.with.HEC$plot.list, function(x){plot(x)})

# show the differences between the MLE and moment estimates (?) from HEC
compare.with.HEC$diff.stat[ grepl('-q$', compare.with.HEC$diff.stat$Measure),] 
##    Year    Measure  value.Q value.HEC          diff     mean        pdiff
## 1  0.01 Q001-avg-q 2.416806     2.438 -0.0211940053 2.427403 0.0087311441
## 2  0.01 Q003-avg-q 2.581846     2.544  0.0378463612 2.562923 0.0147668730
## 3  0.01 Q007-avg-q 2.675951     2.612  0.0639508081 2.643975 0.0241873688
## 4  0.01 Q015-avg-q 2.838811     2.738  0.1008114576 2.788406 0.0361537980
## 5  0.05 Q001-avg-q 2.988016     2.976  0.0120155866 2.982008 0.0040293612
## 6  0.05 Q003-avg-q 3.145622     3.108  0.0376215530 3.126811 0.0120319251
## 7  0.05 Q007-avg-q 3.285016     3.228  0.0570163190 3.256508 0.0175084220
## 8  0.05 Q015-avg-q 3.450660     3.378  0.0726602881 3.414330 0.0212809790
## 9  0.10 Q001-avg-q 3.308550     3.290  0.0185499803 3.299275 0.0056224414
## 10 0.10 Q003-avg-q 3.462336     3.433  0.0293360484 3.447668 0.0085089539
## 11 0.10 Q007-avg-q 3.620229     3.579  0.0412285130 3.599614 0.0114535920
## 12 0.10 Q015-avg-q 3.787091     3.740  0.0470913506 3.763546 0.0125124961
## 13 0.20 Q001-avg-q 3.711903     3.693  0.0189030938 3.702452 0.0051055614
## 14 0.20 Q003-avg-q 3.861207     3.847  0.0142065088 3.854103 0.0036860737
## 15 0.20 Q007-avg-q 4.036056     4.020  0.0160556821 4.028028 0.0039859909
## 16 0.20 Q015-avg-q 4.204141     4.192  0.0121414722 4.198071 0.0028921552
## 17 0.50 Q001-avg-q 4.530930     4.530  0.0009301536 4.530465 0.0002053108
## 18 0.50 Q003-avg-q 4.672124     4.689 -0.0168756631 4.680562 0.0036054778
## 19 0.50 Q007-avg-q 4.862202     4.893 -0.0307980917 4.877601 0.0063141885
## 20 0.50 Q015-avg-q 5.031834     5.074 -0.0421657228 5.052917 0.0083448277
## 21 0.80 Q001-avg-q 5.413735     5.437 -0.0232648779 5.425368 0.0042881662
## 22 0.80 Q003-avg-q 5.547480     5.576 -0.0285196961 5.561740 0.0051278369
## 23 0.80 Q007-avg-q 5.729219     5.766 -0.0367809254 5.747610 0.0063993431
## 24 0.80 Q015-avg-q 5.899323     5.938 -0.0386773799 5.918661 0.0065348189
## 25 0.90 Q001-avg-q 5.901163     5.932 -0.0308372016 5.916581 0.0052119965
## 26 0.90 Q003-avg-q 6.031295     6.047 -0.0157054324 6.039147 0.0026006043
## 27 0.90 Q007-avg-q 6.198934     6.207 -0.0080657874 6.202967 0.0013003112
## 28 0.90 Q015-avg-q 6.368843     6.367  0.0018431866 6.367922 0.0002894487
## 29 0.95 Q001-avg-q 6.317356     6.348 -0.0306438498 6.332678 0.0048390033
## 30 0.95 Q003-avg-q 6.444656     6.435  0.0096555098 6.439828 0.0014993429
## 31 0.95 Q007-avg-q 6.595474     6.559  0.0364736611 6.577237 0.0055454383
## 32 0.95 Q015-avg-q 6.764989     6.704  0.0609888210 6.734494 0.0090561841
## 33 0.99 Q001-avg-q 7.131889     7.140 -0.0081113549 7.135944 0.0011366898
## 34 0.99 Q003-avg-q 7.254258     7.154  0.1002581369 7.204129 0.0139167602
## 35 0.99 Q007-avg-q 7.360602     7.174  0.1866024313 7.267301 0.0256769915
## 36 0.99 Q015-avg-q 7.528800     7.282  0.2467995701 7.405400 0.0333269745